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SIMULATIONS OF MHD INSTABILITIES IN INTRACLUSTER MEDIUM INCLUDING ANISOTROPIC 

THERMAL CONDUCTION 
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ABSTRACT 

We perform a suite of simulations of cooling cores in clusters of galaxies in order to investigate 
the effect of the recently discovered heat flux buoyancy instability (HBI) on the evolution of cores. 
Our models follow the 3-dimensional magnetohydrodynamics (MHD) of cooling cluster cores and 
capture the effects of anisotropic heat conduction along the lines of magnetic field, but do not account 
for the cosmological setting of clusters or the presence of AGN. Our model clusters can be divided 
into three groups according to their final thermodynamical state: catastrophically collapsing cores, 
isothermal cores, and an intermediate group whose final state is determined by the initial configuration 
of magnetic field. Modeled cores that are reminiscent of real cluster cores show evolution towards 
thermal collapse on a time scale which is prolonged by a factor of ~ 2 — 10 compared with the zero- 
conduction cases. The principal effect of the HBI is to re-orient field lines to be perpendicular to 
the temperature gradient. Once the field has been wrapped up onto spherical surfaces surrounding 
the core, the core is insulated from further conductive heating (with the effective thermal conduction 
suppressed to less than 10"^ of the Spitzer value) and proceeds to collapse. We speculate that, in real 
clusters, the central AGN and possibly mergers play the role of "stirrers," periodically disrupting the 
azimuthal field structure and allowing thermal conduction to sporadically heat the core. 
Subject headings: conduction - convection ~ galaxies: clusters: general - instabilities - MHD - plasmas 



1. INTRODUCTION 

The temperature structure of the intracluster medium 
(ICM) in central regions of galaxy clusters is bimodal. 
The non-cooling core clusters have isothermal ICM cores 
with low densities and hence long radiative cooling times. 
On the other hand, the central regions of the ICM in cool- 
ing core clusters have radiative cooling times that can be 
as short as 10® years. Within the cooling radius (the lo- 
cation where the radiative cooling time is comparable to 
the Hubble time and hence the age of the cluster), the 
temperature decreases with decreasing radius (Peterson 
& Fabian 2006). Given the short cooling time, it is pro- 
foundly puzzling that the cores of these clusters have not 
undergone a cooling catastrophe. The current paradigm 
is that the central active galactic nucleus (AGN) heats 
the ICM; this remains one of the most direct arguments 
for "AGN feedback". However, early work (Binney & 
Cowie 1981) highlighted the possible role that thermal 
conduction may play in these clusters. 

Our understanding of the action of thermal conduction 
in atmospheres such as the ICM is undergoing a revolu- 
tion. Because the ICM plasma is very dilute, thermal 
conduction will act in a very anisotropic manner, oc- 
curring essentially unchecked along magnetic field lines 
but being very strongly suppressed perpendicular to field 
lines. This elementary fact has profound implications for 
the dynamics and structure of the ICM ( or, indeed, an y 
dilute plasma atmosphere) . As shown by Balbus ( 2000 1 , 



the anisotropy of conduction fundamentally alters the 
Schwarzschild criterion for convection — rather than re- 
quiring an inverted entropy gradient, convection in such 
an magnctohydrodynamic (MHD) atmosphere will oc- 
cur whenever the temperature gradient is inv erted, due 



to the magneto-thermal instability, (MTI; see Parrish & 
|Stone| [2005 2 0071 for the fi r st num erical studies of this 



instability). [Parrish et al. (12008) studied the effect of 



the MTI in the outer regions of clusters, where tempera- 
ture decreases with radius, and found that the tempera- 
ture profile of the ICM can be substantially modified on 
timescales of several billion years. The instability drives 
field lines to become preferentially radial leading to con- 
duction at a high fraction of the Spitzer conductivity. 

The ICM cores of cooling core galaxy clusters will be 
stable to the MTI (since the temperat ure in suc h cores 
is increasing with radius). However, Quataert (20081 
has discovered a related instability (the heat jiux buoy- 
ancy instability; HBI) which acts when the tempera- 
ture is increasing with radius. Based on local simula- 
tions of the HBI in 3D stratified atmosphere, |Parrish| 



fc Quataert] ( |2008| ) found the HBI induces MHD turbu- 
lence and can somewhat amplify the magnetic field in the 
plasma. They also discovered that, in the plane-parallel 
geometry that characterizes all local simulations, the in- 
stability saturates when the lines of magnetic field are 
re-oriented in such a way as to suppress the heat trans- 
port across the temperature gradient. 
Guided by the results from the local simulations, |Bal-| 
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bus fc Reynolds (20081 suggested that MHD turbulence 
driven by the HBI may be important in regulating the 
conduction of heat into ICM cores and could mediate 
the stabilization of the cooling cores. They hypothesized 
that the presence of radiative cooling and the spherical 
(as opposed to planar) geometry would prevent field line 
re-orientation from completely insulating the core from 
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the conductive heat flux. They also highlighted the fact 
that HBI driven turbulence would create a convective 
heat flux that removes heat from the cool core (i.e., it 
acts as a cooling term in the energy equation) . 

In this paper, we present global models of cooling core 
clusters in which we explore the role of heat conduc- 
tion and the HBI on the evolution of these cores. The 
results of our simulations suggest that HBI alone can- 
not regulate and stabilize a cooling core. We have fol- 
lowed the non-linear evolution of the HBI in the inner 
~ few X 100 kpc in clusters, and found that it is ubiq- 
uitous and rearranges the lines of magnetic field in such 
a way that they are wrapped around the core. This re- 
sults in dramatic suppression of the heat conduction be- 
low the Spitzer value. Consequently, within context of 
these simple models (which do not include AGN or real- 
istic cluster/dark matter dynamics) heat conduction can 
significantly delay but cannot prevent the catastrophic 
core collapse in real clusters. In § 2 we briefiy review 
the equations and time scales governing the problem of 
heat conduction in cluster cores as well as the numerical 
setup used in simulations. In § 3 we present results from 
a parameter space study of cooling core clusters and then 
focus on a specific model of a core resembling that in the 
Perseus cluster. We discuss our results, approximations, 
and observational consequences in § 4, and present our 
conclusions in § 5. 

2. SIMULATIONS 



Using the 3-dimensional MHD code Athena (Stone 



et al.||2008 l, we have carried out simulations oi ther 
mally conducting ICM cores which incorporate the ef- 
fects of anisotropic thermal conduction. We run two 
classes of simulations. Firstly, we conduct a suite of 
simulations aimed at mapping out the behavior of clus- 
ters as a function of position in the 2-dimensional pa- 
rameter space (tcooi/idyn, icond/idyn), whcrc the cooling 



timescale tr 



conduction timescale tr 



cal timescale tdyn are defined in Section 2.3 



and dynami- 
These sim- 



ulations describe "theorist clusters" in the sense that no 
particular physical parameters (e.g., density, tempera- 
ture and timescales) are implied, with only dimensionless 
ratios being relevant. We shall refer to these as our "pa- 
rameter space survey simulations" . Secondly, we shall 
conduct simulations that are specifically tailored to de- 
scribe physically-realistic clusters. We shall refer to these 
as our "physical clusters/simulations" . 

We must note one pecuharity of these (or indeed any) 
simulations of "idealized" galaxy clusters. In remov- 
ing our model clusters from their cosmological setting, 
the computational complexity of the problem is reduced 
enormously. Indeed, it is currently infeasible to per- 
form a cosmological simulation which resolves the small 
scale physics relevant to our study. However, it is im- 
portant to realize that galaxy clusters are dynamically 
young objects and, through the accretion of subclusters 
and groups, are still in the process of forming. Any ab 
initio model for the thermodynamic state of the ICM 
must acknowledge the cosmological setting. Neglecting 
the cosmological setting has two consequences for our 
work. Firstly, real systems will possess dynamics related 
to merging substructures that is not captured in our 
treatment. The effects of neglecting this phenomenon 
will be discussed in Section 4. Secondly, there is no 



well-defined choice of "physical initial conditions" for our 
model clusters. We must be content with forming well de- 
fined initial conditions that describe gross aspects of the 
systems under consideration. We describe our choice of 
initial conditions below. The fact that our results agree 
well with those of Parrish et al. (2009) who employ a 
rather different choice of initial condition suggests a ro- 
bustness to these choices. 



2.1. Equations 
The fundamental equations of the analysis are 

dp 
at 

dv 



+ V • (pv) = 



(1) 



, (V X B) X B „ 
+ /o(v • V)v - ^ An ^P + Z'g' (2) 



OB 

Ik 

de 
di 



V X (v X B), 



(3) 



+ V-(ev)==-pV-v-V-Q-n2A(r), (4) 



where p is the mass density, v is the fiuid velocity, B is 
the magnetic field vector, g is the gravitational acceler- 
ation, p is the gas pressure, e is internal energy density, 
Q is the heat fiux, and A(T) is a cooling function, is 
the electron number density, and T is the temperature. 
We adopt an equation of state p = (7 — l)e = pT/pirip, 
with 7 = 5/3, adequate for monoatomic gas. We have 
ignored viscous terms in the equation of motion in this 
work; these will be considered in a future study. 

In the absence of magnetic fields, we take the elec- 
tron thermal c onductivity to be given by its Spitzer value 
( |Spitzer|[l962l ): 



1.84 X 10-5T5/2 



erg s ^ cm ^ K ^ , 



(5) 



where the Coulomb logarithm In A is 30 — 40 for condi- 
tions in a cluster cooling fiow. The local conductive heat 
fiux will then be given by Q = — xVT. If such efficient 
conduction operated unimpeded, radiative loses within 
massive ICM cores would be more than balanced by the 
conduction of heat from the outer portions of the ICM 
in all but the lowest mass clusters. This would elimi- 
nate the possibility of a cooling flow or , indeed, any sig- 
niflcant departure from isothermality ([Binney fc Cowie 



T98T| [Fabian et ar|[2002l [Voigt et al.||Wi| ). Ho- 



magnetic fields of any plausible astrophysical strength 
will strongly suppress transport of electrons across their 
line of force. Thermal conduction is expected to be ef- 
ficient along field lines, and suppressed perpendicular to 
the field lines. The form of the heat fiux under these 
conditions is, 



Q = -xb(b.vr) 



(6) 



where b is a unit vector in the direction of the mag- 
netic field. Assuming a spherical temperature gradi- 
ent, the radial heat flux can be expressed as Q • f = 
— x(b • f)^ dT/dr. For convenience we will define a ther- 
mal diffusivity, k = xT/p = Kanisoino/n){T/To)^/'^, 
where Kaniso has dimensions of a diffusion coefficient 
(cm^ s^^), and uq and Tq are fiducial values for the elec- 
tron number density and temperature, respectively. 
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Fig. 1. — Examples of the geometry of magnetic field lines immediately after the initial radiative relaxation phase in which the cool 
core is imprinted into the initial conditions. Left: Large scale, loosely tangled magnetic field lines used as a part of initial conditions in 
T-models. Right: Tangled azimuthal field lines with spatially variable {j> and 6 components used in A-models. Fields lines in AR-models 
appear very similar except for the addition of the split monopole term at the center. Both panels show the slice through the xy-plane of 
a cluster center. Arrows mark the direction of magnetic field vectors and colors indicate field strength. The orientation of vectors in the 
corners of figure representing A-model is the visualization artifact that arises for very weak fields. 

less otherwise stated the nominal numerical resolution 
used is 100^. 

All of our simulations have an initial density dis- 
tribution that is described by a /3-model, p{r) = 

po (l + (r/ro)^) ' , where = + y'^ + . The 
temperature distribution is initially isothermal with tem- 
perature Tq. We choose to work in units where pq — 1, 
ro — 0.1 and Tq = 1. The underlying gravitational po- 
tential is assumed to be dominated by the dark matter, 
and is static throughout the simulations with the form 



Anisotropic thermal conduction is implemented in the 
Athen a code via a split ope rator approach with subcy- 
cling ( Parrish fc Stone|[2005 ) . Subcychng ensures stabil- 
ity of the integration and is used whenever the Courant 
time step for conduction falls below that used for hy- 
drodynamic equations. The thermal conduction time 
step is evaluated as Stcond = 0-5 (Sx)'^ /umax, where 
i^max = i^aniso for Tq = 1 and 6x is the size of a resolution 
element. The thermal conduction term is implemented 
using the method of monotonize d central difference lim- 
iters (Sharma fc Ham mett|2007 l which prevents unphys- 
ical heat How from colder to hotter regions that may arise 
as a numerical artifact in a variety of algorithms. 

In this work, we adopt two forms for the optically-thin 
radiative cooling function; both are introduced into the 
Athena code using operator splitting as explicit source 
terms. The majority of our parameter space survey sim- 
ulations employ a very simple "bremsstrahlung" form 
for the cooling function, n\K = a'n? T°-^. We control 
the strength of radiative cooling and hence the cool- 
ing timescale via the parameter a. On the other hand, 
our p hysical cluster simulations use the [Tozzi fc Norman] 
(2001| approximation for the cooling function that incor- 
porates the effects of free-free and line cooling: 
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In 



1+1^ 



(8) 



nlA= (Ci(fcsr)-i-^ + C2(A:bT)"-5 + C3) n^n^, (7) 

where (ksT) is in units of keV and is the ion num- 
ber density. For mean metallicity of Z = 0.3Zq, Ci = 
8.6 X 10-3, C2 = 5.8 X 10-2, and C3 = 6.4 x 10"^ and the 
units of A are lO-^^ergcm^^ s-^. Mean molecular weight 
corresponding to this metallicity in case of a near com- 
plete ionization is p — 0.59. 

2.2. Numerical Setup and Initial Conditions 

The simulations were performed in a Cartesian coordi- 
nate system (x, y, z) with a cubic spatial domain defined 
by a; = ±L/2, y = ±L/2, z = ±L/2, where L = I. Un- 



where is the speed of sound corresponding to temper- 
ature Tq. This form corresponds to hydrostatic equilib- 
rium in the initial isothermal ICM density distribution. 

The simplicity of an isothermal hydrostatic atmosphere 
makes it a compelling choice of initial condition. How- 
ever, the astrophysical systems of interest have ICM cores 
that display a positive temperature gradient. Further- 
more, the interesting dynamics is driven by tempera- 
ture gradients. Thus we employ a pre-cursor simulation 
where, starting from the isothermal state, we create a 
cool core in the atmosphere by allowing radiative cool- 
ing in the absence of conduction. This continues until the 
core reaches a temperature of Tc ~ 0.3 — 0.4To, similar to 
the range observed in cores. The density and magnetic 
field are allowed to evolve self-consistently during this 
initial cooling phase. We then use this cool-core state 
as the starting point for the full simulations (i.e., those 
including the anisotropic conduction and associated dy- 
namics). 

Note that cluster cores initialized in this way are not 
in thermodynamic equilibrium and that anisotropic con- 
duction is not matched to exactly balance radiative cool- 
ing at the start of the full simulation (for alternative 
approach starting from thermodynamic equilibrium see 
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TABLE 1 
Value of Kaniso in models. 



Model 


B-field 


Bl'' 


B2 




B3 


B4 




structure 


(a=0) 


(o=0.01) 


(a Si 0.1) 


(a=0.3) 


(«=1) 


1 


T 








0.1 




2 


T 


0.025 


0.025 




1 


1 


3 


T 


0.05 


0.05 




10 


10 


4 


A 


0.01 










5 


A 


0.025 


0.025 


0.0 


0.1 




6 


A 


0.1 


0.05 


0.1 


1 


1 


7 


A 


1 




1 


10 


10 


8 


AR 


0.01 


0.01 


0.01 






9 


AR 


0.025 


0.025 


0.025 


0.1 




10 


AR 






0.1 


1 


1 


11 


AR 






1 


10 


10 



^ Bl— 4: Models with bremsstrahlung cooling function, n^A = an^T^'^, 
and fi = 1. 

^ TN: Models with Tozzi-Norman cooling function for metallicity Z = 
0.3 Zq and = 0.59. 



[Parrish et al. 2009). While realistic cluster cores are 



most likely not m thermodynamic equilibrium either, 
they experience more gradual thermal histories compared 
to that created in our simulations. We will discuss the 
effect of this feature later, along with the description of 
results. 

The initial magnetic field strength is typically chosen 
so that the average value of the plasma parameter (i.e. 
the ratio of thermal pressure to magnetic pressure) is 
P — Snp/B'^ ^ 10^ — 10^. We explore several differ- 
ent scenarios for the initial structure of magnetic fields. 
Firstly, we define field line geometry given by the form. 



— ^ cos(7r y) cos(7r z) , 
v3 

By = cos(7r x) cosfTT z), 

Bz — —2: cos(7r x) cos(7r y), 
V3 



(9) 
(10) 

(11) 



This initial geometry is self-consistently evolved within 
an isothermal core until field lines assume loosely tangled 
geometry while retaining the large scale character given 
by the initial form. This field configuration is then used 
as a part of the initial conditions. These models are rep- 
resentative of an idealized scenario in which anisotropic 
conduction initially operates efficiently by transporting 
heat from the outer, hotter regions of a cluster to the 
cooler core center. Hereafter, we denote this class of 
models by "T" . In this model, the constant Bq is the ini- 
tial amplitude of the magnetic field. Secondly, we exam- 
ine models with tangled purely azimuthal field lines, i.e., 
a field structure that is initially wrapped onto surfaces of 
constant r. In this scenario we investigate the evolution 
of cool cores when the initial field geometry is set up to 
inhibit heat conduction towards the core. We use these 
models to test the hypothesis that spherical collapse and 
MHD turbulence driven by the heat fiux instability can 
regulate field-line insulation and drive a rev erse convec- 
tive thermal flux ( Balbus & Reynolds 2008 1 . Defining a 
spherical polar coordinate system (r, 0) where 9 = 
is aligned along the z-axis, we employ an azimuthal field 
structure defined by. 



B^ = 2Bo(l -I- sin(27rr/r2)) sin(36l) - 

Bo{l + sin(27rr/ri)) sin(2(?!)) sin(26l), (13) 

Br = 0. (14) 
We refer to these as "A" models. Here, Bq = Snp/f], 
ri = 0.1 and r2 = 0.087 are coherence lengths defin- 
ing characteristic scales on which magnetic field vector 
changes direction. Finally, we examine a modification of 
the tangled azimuthal field model in which we include a 
radial component in the form of a split monopole, 

Bro . f 1 if z > 



Bj. = sign 



sign = 



' \ otherwise 



(15) 



= 2Bo(l + sin(27rr/ri)) sin6'cos(2(?!i). 



(12) 



We shall refer to these as our "AR" models. The constant 
e << ro is a small number chosen to avoid a singularity 
at the origin. We use this simple setup to study cases 
in which magnetic field at the very center of a cool core 
may be tangled and enhanced in the process of spherical 
collapse. We test to which extent such fields provide a 
local magnetic support as well as channel the heat within 
a small volume of the core. All initial magnetic field 
geometries are self-consistently evolved during the initial 
precursor simulation in which the cool core is imprinted 
into the initial condition. For illustration, slices through 
the "T" and "A" field geometries (at a time just after 
the initial radiative relaxation of the model cluster) are 
shown in Fig. [T] 

We complete the specification of our simulations by de- 
scribing the boundary conditions. The velocity bound- 
ary conditions at the interface of the active and ghost 
zones are zero-gradient outfiow, and the pressure and 
density of the gas are extrapolated from the active zone 
into the ghost zone so as to maintain hydrostatic equilib- 
rium in the ghost zones. This choice of boundary condi- 
tions ensures stability and prevents spurious sound waves 
from developing at the boundaries. The temperature is 
fixed to the virial temperature Tq on a spherical shell 
at a radius r = L/2. Thus, while the dynamics are fol- 
lowed in the full cubical domain, only the interior of the 
r = L/2 sphere is physically interesting. Early exper- 
imentation with applying temperature boundary condi- 
tions at the edge of the cubic domain revealed behavior 
that appeared to depend upon the overall orientation of 
our Cartesian domain. Thus, we chose to impose ther- 
mal boundary conditions on a spherical surface in order 
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Fig. 2. — Illustration of the parameter space of clusters explored 
in this study shown in terms of cooling and heat conduction time, 
normalized to the dynamical time. Also marked are corresponding 
values of Kaniso and a (for a list of values used in simulations see 
Table [ij. Filled circles, diamonds, and squares represent T, A, 
and AH-models, respectively. A combination of symbols signifies 
that different magnetic field configurations were simulated at the 
same value of a and Kaniso- Letter "P" marks the location of the 
Perseus-like model. The investigated parameter space includes but 
it is not limited to that of realistic clusters which occupy the range 

1 £ <cool/*dyn < 10 and few X 0.1 < tcond/^dyn < 10. 

to emulate hot intercluster plasma surromiding the cool 
core. 

2.3. Characteristic Time Scales 

For our parameter space survey, we catalogue our simu- 
lated clusters according to their characteristic dynamical, 
conduction, and cooling timescales. These are defined (in 
code units) by. 



^dyn — 

^cond — 
^cool — 



R _ 0.39 

To 0.036 



K AT ftaniso M 
e U 

^ ^0-27- 

nth. a 



(16) 

(17) 
(18) 



where ^dyn is calculated for the cooling core radius R — 
L/2 — 0.5 and the average speed of sound in the region 
with virial temperature Tq = 1, 7 = 5/3, and /i = 1 
or 0.59, as specified for a given run. When calculating 
icooi we evaluated the shortest time scale corresponding 
to the core center where typically Tc « 0.3 and pc ~ 
3 at the beginning of a simulation. In eq. (17 1, k was 



evaluated for an average density in the core region (p) — 
0.1 and virial temperature.^ According to the theory of 
linear growth applied to the weak magnetic field regime, 
the plasma is unstable to HBI on the local dynamical 
time whenever the vector of magnetic field is aligned with 
the temperature gradient and unbridled conduction is 
allowed in the direction of gravity. In the limit when 
lines of magnetic field are preferentially oriented across 

Note that the expression commonly used in literature to es- 
timate tcond omits factor To/ AT and thus, refers to the shortest 
timescale for heat conduction. 



the temperature gradient, this condition is modified by 
a factor (b • f ) in such way that 

-1/2 



HBI 



dlnTd^ 
dr dr 



(19) 



( Quataertj|2008[ ). For the values of initial magnetic field 
used in simulations the Alfven time scale is typically 
much longer with respect to the above time scales and 
thus, Alfven waves are not expected to play a significant 
role during the phase characterized by the HBI instabil- 
ity and radiative cooling. 

3. RESULTS 

We have performed an extensive suite of simulations. 
Our primary set of 36 simulations (our parameter space 
survey simulations) represent a systematic exploration 
of the dependence of the cluster core evolution on ra- 
tios of the characteristic timescales (i.e., tcooi/^dyn and 
icond/^dyn), and the magnetic field geometry. The cool- 
ing and conduction timescales are controlled via the pa- 
rameters a and Kaniso- Table 1 defines this set of mod- 
els. We explore three field geometries (the T, A, and 
AR geometries of Section 2.2) and five cooling laws (Bl- 
4-l-TN; this includes the zero-cooling case a = 0). For 
each choice of field geometry and cooling law, 2-4 sim- 
ulations are performed with different degrees of thermal 
conduction. Thus, we can map out the behavior of clus- 
ters as a function of their position on the 2-d parame- 
ter space (tcooi/^dyn, ^cond/idyn) and, furthermore, assess 
whether the initial field geometry can infiuence/change 
the eventual state of the cluster for a given position on 
the (^cooi/^dyn, icond/idyn)-plane. We refer to individual 
models from Table [T] according to their alphanumeric tag 
consisting of the letter abbreviation corresponding to the 
magnetic field structure (T, A, or AR), a model num- 
ber (1—11), and the descriptor of the cooling function 
(Bl— B4 or TN). According to this convention a model 
with purely azimuthal initial magnetic field geometry, 
f^aniso = 0.025, and no radiative cooling is marked as 
A5B1. Note that this study includes but it is not limited 
to the part of the parameter space occupied by realistic 
clusters. 

We expect the HBI to be active; hence, one could read- 
ily envisage a scenario in which field line reorientation 
(by the HBI) completely insulates the central region from 
the conductive heat fiux, which would prevent complete 
equilibration of core. However, as we shall see, while 
field line reorientation does occur, it never completely 
insulates the core. Figure |2] shows the models on the 
2-d parameter space (^cooi/^dyn, ^cond/^dyn)- The results 
of our simulations clearly divide this parameter space 
into three regions depending upon the ratio tcooi/^cond- 
If icooi/^cond ^ 10^, coohng appears to be too weak to 
sustain any significant temperature gradient against the 
action of conduction, and all simulated cluster cores be- 
come approximately isothermal. On the other hand, if 
icooi/icond ^ 25, our simulated clusters always undergo 
a cooling catastrophe although it can be appreciably de- 
layed by the action of HBI-rcgulated conduction. In the 
intermediate regime (25 < tcooi/^cond ^ 10^), clusters 
can evolve to either an isothermal state or undergo a 
cooling catastrophe depending upon the initial magnetic 
field configuration. These are the only outcomes; none of 
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TABLE 2 

Time scale (in tdyn) for core evolution towards collapse (C) or 

ISOTHERMAL STATE (I). 



Model 


B-field 


Bl 


B2 


TN 


B3 


B4 




structure 


(a=0) 


(a=0.01) 


{a X 0.1) 


(a=0.3) 


(a=l) 


1 


T 








12 (C) 




2 


T 


70 (I) 


> 1000 (C) 




80 (C) 


10 (C) 


3 


T 


60 (I) 


80 (I) 




<3 (I) 


<3(I) 


4 


A 


> 500 (I) 










5 


A 


500 (I) 


380 (C) 


30 (C) 


8(C) 




6 


A 


240 (I) 


510 (C) 


160 (C) 


18 (C) 


5(C) 


7 


A 


50 (I) 




> 700 (I)^ 


20 (I) 


13 (C) 


8 


AR 


> 500 (I) 


420 (C) 


90 (C) 






9 


AR 


500 (I) 


560 (C) 


110 (C) 


8(C) 




10 


AR 






230 (C) 


20 (C) 


5 (C) 


11 


AR 






> 700 (V)^ 
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A7TN: Evolution of the core that initially headed towards collapse was reversed 
to evolution towards isothermal state. 

^ ARllTN: Final state for this model appears undetermined (U) at the end of 
the simulation. 



our simulated clusters achieved a non-isothermal quasi- 
steady state. 

In the rest of this section we shall describe these results 
in detail and, in particular, examine the role of the HBI 
in influencing the thermal evolution of the cluster core. 

3.1. Clusters with an isothermal final state 

As an extreme case, we first describe the evolution of 
cool cores in models without radiative cooling, a — 
(also referred to as "Bl" models). The fact that these 
clusters evolve to an isothermal state is not as trivial 
as it first appears. Recall that our full-up simulations 
start with a cool core (and hence a positive temperature 
gradient) already imprinted on the ICM atmosphere. We 
expect the HBI to be active and, hence, one could readily 
envisage a situation in which field-line re-orientation by 
the HBI completely insulates the core from the conduc- 
tive heat flux and hence prevents complete equilibration 
of the core. As we will see, field-line re-orientation does 
indeed occur but never completely insulates the core. 

Table [2] (Bl column) shows the time taken for each 
of the non-cooling models to achieve an approximately 
isothermal state. Cores with loosely tangled field lines, 
captured in T-models, evolve to isothermal state on sig- 
nificantly shorter time scales with respect to the cores in 
models A and AR, where the field line geometry is pref- 
erentially azimuthal. Figure |3] shows the evolution of the 
core temperature and mean angle between magnetic field 
lines and the radius unit vector, f, in three models repre- 
senting different magnetic field geometries, T2B1, A5B1, 
and AR9B1, all at the same value of Kaniso = 0.025. 
The core temperature was measured at the center of 
the core, whereas the angle {9b) is calculated as a vol- 
ume average of cos^^(5 • f), and it indicates the align- 
ment of the field lines with the global temperature gra- 
dient^. As expected, the action of the HBI is most ap- 
parent in the tangled field model (T2B1) where we see a 
rapid reorientation of the field lines from (^b) ~ 57° to 
{6b) ~ 71°. The corresponding effective thermal conduc- 

^ The region of the cluster core outside of 0.3 L is not taken 
into account in the calculation of (6b), in order to focus on the 
field structure closer to the core center and avoid effects of the 
boundaries. 



tivity is (b • f )^ w 0.25 of the Spitzer value but decreases 
as the field lines arc re-oriented. Once the core has equi- 
librated, the HBI is no longer driven and the field lines 
become more disoriented again leading to a final value 
between {9b) ~ 60° — 70°. A similar final field orienta- 
tion results in the A- and AR-models (i.e., the azimuthal 
initial field configurations). In these cases, the core takes 
substantially longer to equilibrate, however, with the ini- 
tial heat fiux corresponding to only ^ 0.02 of the Spitzer 
value. In comparison, the time scale for unbridled con- 
duction in this set of models is icond ~ 3.7tdyn- 

Given the same initial geometry of magnetic field lines, 
the cores with higher values of Kaniso evolve towards 
isothermality on shorter time scales. The behavior of 
models shown in Figure |4] and Fi gure \5\ confirms this 
general expectation, summarized in eq. p7|. Figure [i] 
shows the evolution of the core temperature and orien- 
tation of the magnetic field lines in T-models, T2B1 and 
T3B1. Similarly, Figure [s] shows the evolution for four 
A-models, A4B1, A5B1, A6B1, and A7B1. 

Interestingly, these non-cooling cores remain rather 
kinematically quiescent throughout the entire evolution, 
and magnetic field does not play a significant role in 
subsequent dynamical evolution of the core. Figure |6] 
shows that the early evolution of the magnetic and ki- 
netic energy density exhibit different behaviors in T- and 
A-models. In the T-models, the volume averaged mag- 
netic and kinetic energy increase over time due to the 
action of the HBI; it appears that a weak and short lived 
dynamo may be established. The A- and AR-models, on 
the other hand, show a monotonic decrease in the kinetic 
and magnetic energy densities as the non-equilibrium az- 
imuthal magnetic fields relax and numerically reconnect. 
In both panels the components of energy density are nor- 
malized to the internal energy density of the gas and, 
hence, the kinetic and magnetic energy represent a small 
fraction (~ 10~*) of the internal energy. This implies 
that cores are kinematically quiescent and that even af- 
ter it has been enhanced (T-models), magnetic field is 
not sufficiently strong to be dynamically important. 

In fact, these qualitative results are borne out even 
when there is cooling present provided it is weak and 
much slower than the heat conduction in the sense that 



7 



1.0 - 



i: 0.8 

^ 0.6 

s 

oj 0.4 

o 
u 

0.2 - 




0.0 



100 200 300 400 
Time (t,yj 



500 



600 



90 



80 



70 



60 



50 



100 200 300 400 
Time (tjyj 



500 



600 



Fig. 3. — Left: Evolution of the core temperature in T2B1 (solid), A5B1 (dashed), and AR9B1 (dash-dot) models, respectively. Right: 
Evolution of the mean angle, {9g), between the lines of magnetic field and the radius unit vector (same models), indicating the alignment 
of the field lines with the global temperature gradient. 
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Fig. 4. — Left: Evolution of the core temperature in T2B1 (solid) and T3B1 (dashed) models, respectively. Right: Evolution of the ' 
for the same two models. 
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icooi/idyn > lO'^. Examples of such behavior are models 
T3B3, A7B3, and AR11B3, all of which evolve towards 
isothermal final state. 

3.2. Clusters which undergo catastrophic collapse 

Real cooling core galaxy clusters do not have isother- 
mal ICM atmospheres and hence, while this avoids the 
cooling catastrophe, it must not be considered a physical 
outcome. The more interesting case are those clusters in 
which cooling is too effective to allow an isothermal state 
to be achieved. As already mentioned, all of our model 
clusters with tcooi/^dyn < 25 undergo an eventual cooling 
catastrophe; Table |2] lists the time taken for the model 
clusters to either undergo catastrophic collapse (C) or 
equilibrate to an isothermal state (I). We now discuss 
the collapsing clusters in more detail. 

To illustrate many of the dynamical aspects of these 
collapsing clusters, we discuss run T2B2 in detail. Fig- 
ure n\ (solid line) shows the evolution of the core tem- 
perature together with the average field line orienta- 
tion for run T2B2. At t = 0, the core is not in a 
state of thermal balance, with conductive heating over- 
whelming the radiative cooling. The core is rapidly 



heated from T w OATq almost to a state of isother- 
mality (T w 0.85To). The HBI is clearly playing a 
role during this early heating event, as revealed by the 
rapid re-orientation of the field geometry towards an az- 
imuthal configuration (Fig. [8|. After this early heating 
event, an approximate balance between conductive heat- 
ing and radiative cooling is maintained for a duration of 
~ 300idyn- Eventually, azimuthal field line wrapping re- 
duces the conductive heat fiux to the point where the core 
temperature starts to decrease. This slow collapse will 
end in a cooling catastrophe, in this case at some time 
t > lOOOfdyn- For comparison, the nominal cooling time 
scale for this group of models set by the value of a and 
calculated according to the equation 18 is icooi ~ 69idyn- 



Thus, we see that thermal conduction can dramatically 
increase the time taken for the cluster to undergo the 
cooling catastrophe even when the HBI is inducing az- 
imuthal field-line wrapping within the cluster core. 

As may be expected, cores with initially loosely tan- 
gled fields (T-models) evolve towards collapse on time 
scales longer than cores with initially azimuthal fields 
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Fig. 5. — Left: Evolution of the core temperature in four A-models: A4B1 (solid), A5B1 (dashed), A6Bl(dash-dot) and A7B1 (dash-3-dot), 
respectively. Right: Same four models, evolution of {9b). 
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Fig. 6. — Evolution of the Isinetic (solid) and magnetic (dashed) energy density in two T-models shown in Figure |4] (left) and four 
A-models from Figure [5] (right) without radiative cooling. Both components of energy density are normalized to the instantaneous internal 
energy density. 



(A- and AR-models)^. Furthermore, the AR-model clus- 
ters (which include a split monopole term to the ini- 
tial field configuration and hence possess a significant 
radial field in the core regions) have a collapse that is 
delayed with respect to A-models. Figure IT] compares 
the evolution of core temperature and (0b) in models 
T2B2, A5B2, and AR9B2, characterized by different ini- 
tial magnetic field structure and identical cooling and 
conduction (Kaniso = 0.025 and a — 0.01). In A and AR- 
models, the field orientation starts at (Ob) ~ 80°; this 
permits sufficient conductive heat ffux to raise the core 
temperature to T w 0.5 — 0.6Tq. However, the action 
of the HBI further increases the field line orientation to 
{9b) > 85°, and insulates the core which proceeds to col- 
lapse. The cores in models A5B2 and AR9B2 collapse in 
380 and 560 dynamical times after the beginning of the 
simulation, respectively, compared with the collapse on 
a timescale of > lOOOtdyn for model T2B2. 

® Even in the A-model, the field configuration at the start of the 
full simulation is not perfectly azimuthal due to field distortions 
imprinted during the precursor simulation in which the cool core 
is formed. 



As was the case for the non-cooling clusters, the mod- 
eled cooling cores remain kinematically quiescent and the 
magnetic field never plays a direct dynamical role. The 
evolution of kinetic and magnetic energy density for mod- 
els T2B2, A5B2, and AR9B2 is shown in Figure [9| The 
kinetic and magnetic energy are a small fraction of the 
internal energy of the gas and the core region. Velocities 
remain very subsonic, v < lO^^Cg. Independent of the 
assumed initial geometry of the field lines, HBI does not 
result in significant field amplification. 

The fact that the cores remain (relatively) kinemati- 
cally quiescent results in very little convective heat ffux. 
To see this. Fig 10 shows the evolution of the radial con- 
vective, conductive and advective heat ffow for models 
T2B2 and A5B2 within a sphere of radius r = 0.1 along 
with the total radiative loses and the fiducial value of un- 
bridled conductive heat luminosity at the Spitzer value, 
for comparison. To construct these quantities, we start 
with the radial components of the advective, convective 
and conductive heat fiuxes defined by 
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Fig. 7. — Left: Evolution of core temperature in models with different magnetic field structure, all with Kaniso = 0.025 and a = 0.01: 
T2B2 (solid), A5B2 (dashed), and AR9B2 (dash-dot). The nominal cooling time scale for this group of models sot by value of a is 
tcool ~ 69fdyn- Right: Same models, evolution of magnetic field orientation. 
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Fig. 8. — Final orientation of the magnetic field lines in model 
T2B2, 800 dynamical times after the beginning of the simulation. 
Large scale, loosely tangled magnetic field lines with (6b) = 58° 
are used as a part of the initial conditions in this model, as illus- 
trated in the left panel of Figure ^ Traces of the initial topology 
are still present at late times in the strongest component of the 
field, however, most of the lines are wrapped around the core as 
a consequence of HBI. Final value of (^b) = 85°- Figure shows a 
slice through the xy-plane of a cluster center. 

Qconv = ^ -, fcg ((I'r) k^nbT) -I- (n) {6nbT) + 
7- 1 



(21) 
(22) 



Here, Vr is the radial component of the velocity field. The 
fluctuations {5vr,5n,5T) were deflned as relative to a 
mean calculated on the surface of the sphere. We then in- 
tegrate these components across the surface of the spher- 
ical shell r = 0.1 to obtain relevant luminosities plotted 
in Fig[TOj The convective heat flux acts as a cooling term 
for the core (Balbus & Reynolds 2008) but remains signif- 
icantly smaller than the other heat fluxes. As expected, 
the conductive heat component initially decreases as the 
HBI wraps the lines of magnetic field around the core. 
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Fig. 9. — Evolution of the kinetic (solid) and magnetic (dashed) 
energy density for models T2B2 (ending at 1000tdyn)i A5B2 
(380tdyn)i and AR9B2 (560tdyn). Both components of energy 
density are normalized to instantaneous internal energy density 
and comprise a small fraction of it. 

At later times (t > 600idyn) it starts increasing again as 
a consequence of radial field line stretching in the pro- 
cess of core collapse, which once again "unlocks" heat 
conduction along the temperature gradient. Note that 
trends in evolution of luminosity components in T and 
A-models appear qualitatively similar, however, there is 
about an order of magnitude difference in the magnitude 
of conductive luminosity between the two. This is a con- 
sequence of an early, efficient suppression of anisotropic 
heat conduction by magnetic field lines in model A5B2 
in comparison to T2B2. In both runs the field lines even- 
tually saturate at an angle of {O-q) ~ 85 — 88°, in such 
way suppressing the conductive heat flux to a fraction 
of only (b • f )^ k, of the Spitzer value throughout 

the core volume. Apart from unbridled heat conduction, 
the largest of the remaining four components of luminos- 
ity in both models is the radiative cooling. In terms of 
luminosity magnitudes, radiative cooling is followed by 
advective luminosity which can be interpreted as heat- 
ing of the core in the process of adiabatic compression, 
which is in turn followed by anisotropic conduction and 
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TABLE 3 

Models of Perseus-like cluster. 



Model 


B-fiold 
structure 


Resolution 


'^aniso 
(KSpitz) 


^^collapse 

(Gyr) 
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Fig. 10. — Evolution of luminosity components in models T2B2 (left panel) and A5B2 (right) measured within the sphere of radius of 
0.1 from the center of the cluster. Each component is marked with a different line: anisotropic heat conduction (solid), radiative (dashed), 
convective (dash-dot), advective (dash-3dot), and fiducial value of unbridled heat conduction at the Spitzer value (thick, long-dashed), all 
in arbitrary units. The heat flux at full Spitzer value was divided by a factor of 5 (left panel) and 10 (right) in order to show it alongside 
of other components. Negative fluxes are inflowing and vice versa. 

find that cores with tcooi/^cond < 25 evolve towards col- 
lapse regardless of their initial field geometry (T2, A6, 
and ARIO in case of both B3 and B4 runs but also T2B2 
and A5B2), while those with tcooi/^cond ^ 250 all evolve 
towards isothermal state (T3B3, A7B3, and AR11B3). 

3.3. Models of physical clusters 

In order to offer a more intuitive interpretation of the 
role of HBI for a portion of the parameter space popu- 
lated by real clusters, we present a second group of cal- 
culations where we scale simulation parameters specif- 
ically to match the cooling core of a rich galaxy clus- 
ter. We mark the position of these models with let- 
ter "P" in (tcooi, icond) space in Figure [2j It is worth 
noting that characteristic timescales of real cooling core 
clusters occupy the range 1 < tcooi/idyn ^ 10 and 
few X 0.1 < icond/^dyn ^ 10, which places them in the 
collapse region of Figure |2] 

Specifically, we consider a cluster characterized by a 
virial temperature of kTo = 7 keV, initial core tempera- 
ture of kTc « 3keV, and core radius Rq = 100 kpc; these 
properties are reminiscent of the Perseus cluster. The gas 
number density in the modeled core is uq = 0.03 cm~'^ 
(somewhat lower than in Perseus) and the speed of sound 
is Cs ~ 900 kms^^ in the core center. The core cools ra- 
diatively according to th e Tozzi-Norman cooling function 
( Tozzi fc Norman||2001 1. For this cluster, we ran a set of 
models spanning different magnetic field configurations 
and simulation resolutions, as detailed in Table 3. In this 
group of models we test the evolution of the cores with 
predominantly azimuthal magnetic field geometry, a sce- 
nario that allows us to calculate a lower limit on the life 
time of conducting cores. With the exception of model- 
P5, the (anisotropic) heat conduction is set to the Spitzer 
value (equation p\ and mediated by magnetic field in 
all models. Moder-P5 is a pure radiative collapse model 

(^aniso 

= 0) given for comparison. The (initial) charac- 



convective luminosity. 

As stated at the beginning of this Section, the pa- 
rameter space is occupied by three distinct groups of 
systems which end up in either catastrophic collapse 
(icooi/^cond < 25), isothermal state (icooi/icond > 10^), or 
as border line cases. Border line cases represent the tran- 
sition population between the collapse and isothermal 
group of objects and their final state is determined by the 
structure of magnetic field. The behavior of borderline 
cases is illustrated by comparing models T3B2 and A6B2 
(which differ only in their initial field structure), where 
the former evolves to isothermal state in 80 tdyn and the 
latter collapses after 510tdyn- Another example is pro- 
vided by models T3B4, A7B4, and AR11B4, of which 
only the core in the T-model evolves towards isothermal 
state while those in A- and AR-models collapse. Cores 
in models ATTN and ARllTN are also border line cases 
given that both evolve on very long time scales and that 
they seem to hang on a very edge between collapse and 
isothermality. Interestingly, the core in ATTN initially 
evolves towards collapse but after ~ 300tdyn the tem- 
perature curve reverses towards isothermal. The core 
temperature in ARllTN reminds relatively flat until late 
time in the simulation and its final state undetermined. 
The cause of different outcomes in these two models is 
probably stochastic, driven by a slightly different evo- 
lution of the magnetic field and it underlines the role 
of the field structure in these transition population of 
objects. The values of tcooi/^cond for TN, B2, and B4 
border line cases are 26, 38 and T5, respectively. We also 



teristic time scales for this system are idyn ~ 1-3 x 10 yr 
icond ~ 4.1 X 10'' yr, and icooi ~ 2 Gyr. We assume ini- 
tial values of plasma parameter and the strength of radial 
component of magnetic field to be /3 = 10"^, Br = and 
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Fig. 11. — Left: Evolution of temperature in the core of a cluster for models PI (solid), P2 (dashed), and P5 (dash-dot) models. Right: 
Same for P5 only, starting from an isothermal core at virial temperature of 7 keV. 
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Fig. 12. — Left: Evolution of {6b) in PI (solid) and P2 (dashed) models in the central 60 kpc region. Right: Same, except for P4 model 
(pure radiative collapse). Spherical collapse in absence of HBI leads to smaller {6b) and does not maintain the azimuthal orientation of 
the magnetic field lines. 



Br = 3fJ.G, for A- and AR-models, respectively. 

The collapse of conducting cores is postponed by a 
factor of ^ 2—10 with respect to the time scale for 
pure radiative collapse. The cores in models PI and 
P2 evolve through an initial heating of the core during 
which time the HBI re-orients the field lines in order 
to reduce the conductive heat flux. Figures [TT] and [12] 
show the evolution of the core temperature and (Ob) in 
PI- and P2-models along with the pure radiative col- 
lapse case, for comparison. The core with initial temper- 
ature of '--^ 3 keV, in absence of heat conduction heads 
towards collapse in only tcooi ~ 2 Gyr (model P5). The 
cores in PI and P2 models evolve through the initial 
rise in temperature, have a phase of uniform evolution 
approximately linear with time, and in the final 2 Gyr 
dominated by radiative cooling behave similarly to the 
pure radiative collapse case. As remarked in the previ- 
ous section, the initial rise in the core temperature in 
models with heat conduction is a consequence of non- 
equilibrium initial conditions. The intermediate phase 
of steady, quasi-linear decay of the core temperature is 
specific to conducting cores and is absent in case of pure 



collapse model. The presence of HBI is indicated by the 
high value of {6b) = 87° shown in the left panel of Fig- 
ure [12] contrary to the pure radiative collapse scenario 
P5 (right panel) which does not maintain the azimuthal 
structure of the field. The field line re-orientation ap- 
pears to saturate at {6b) ~ 87° and, from then on, the 
core undergoes a gradual decline in core temperature un- 
til it reaches a temperature of kT ~ 2 — 3 keV. After that 
time, the core temperature decreases more rapidly due to 
the onset of line-cooling within the TN cooling function 
and the cooling catastrophe is reached ~ 2 Gyr later. In- 
deed, this final stage in the conducting core models is 
very similar to the pure radiative collapse case, conduc- 
tion having very little effect. 

If real cluster cores are similar to the cores modeled 
here, they are likely to be observed precisely during the 
long gradual decline phase (although, compared with 
model clusters during this phase, real clusters have a 
larger fractional temperature decrease as one heads into 
the core). In model-Pl, this phase starts from 1 Gyr after 
the beginning of the simulation (at which time the core 
has a temperature of 5 keV) and the core collapses 8 Gyr 
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Fig. 13. — Evolution of mean Itinetic (solid) and magnetic 
(dashed) energy density in Pl(shorter) and P2 (longer) normal- 
ized to mean instantaneous energy density. 

later. In model P2, the life time of the core is extended 
due to the presence of a split monopole, with collapse 
occurring 9.5 Gyr after the start of this decline phase. 
A pure cooling simulation starting from a temperature 
of 5keV undergoes the cooling catastrophe in approxi- 
mately 4 Gyr. Thus, even starting off with rather az- 
imuthal field geometries, thermal conduction delays the 
thermal collapse of the cluster core by at least a factor 
of 2. 

It follows from the findings in previous section that 
the longest lived cooling cores are those with large scale 
magnetic fields (represented by T-models). In order to 
estimate the factor by which core collapse is postponed 
in these models, we consult our grid of models outlined 
in Table [T| and [2] It is possible to infer that T-model 
cores collapse on time scales 2 — 5 times longer with re- 
spect to the cores with predominantly azimuthal fields 
(A-models), as illustrated in examples of model pairs: 
T2B2 - A5B2, T2B3 - A6B3, and T2B4 - A6B4. Hence, 
if T-model cores evolve 2 — 5 times longer than A-model 
cores, and Perseus-like A-model cores evolve on time 
scales ~ 2icooi, than the overall evolution of collapsing 
cores can be prolonge d by a to tal factor of ~ 2 — 10 icooi- 

As found in Sections |3.1fj3.2| these model clusters cores 
are characterized by an absence of strong turbulence or 
a magnetic dynamo effect. The magnetic and kinetic 
energies uniformly decrease in PI and P2 models (Fig- 
ure 
the 



13 1 and are only a small fraction of ~ 10~^ — 10~^ of 
internal energy. Figure 14 shows the distribution of 
velocity magnitude in the xy-plane containing the core 
center. The highest velocities in the core are of the order 
of 30 km s^^ and are visibly associated with the magnetic 
field structure. As the relative strength of magnetic field 
decreases in the process of core collapse, the core cen- 
ter becomes kinematically more quiescent with respect 
to the initial velocity distribution. 

Modeled Perseus-like cores exhibit the strongest heat 
conduction and heat flux buoyancy instability within the 
radius of ~ 20 kpc - a distance comparable to the half- 
depth radius of the core where the temperature gradient 
is steepest. Figure [15] shows the evolution of five lu- 
minosity components in units of ergs^^ calculated for 
a sphere of radius of 20 kpc in model PI. Note that 



during the first ~ 1 Gyr of evolution, anisotropic heat 
conduction overcompensates for the radiative losses and 
that this phase is coincident with the initial rapid rise 
of the core temperature seen in Figure 111 Between 1 



and 7.5 Gyr radiative losses are comparable to but larger 
than the rate of heat conduction and in this period core 
temperature steadily declines. In the remaining period 
of evolution, the rate of radiative cooling becomes > 2 
times larger than the anisotropic heat conduction and 
the core evolves towards collapse in the fashion similar 
to a pure radiative collapse scenario. The rate of energy 
transported by convection is consistently ~ 2 orders of 
magnitude lower than that of the anisotropic conduction 
and it acts as a cooling term throughout the simulation. 
Also included in Figure 15 is the fiducial rate of unbri- 



dled heat conduction at the Spitzer value (thick, long 
dash line), plotted for comparison - because this compo- 
nent of luminosity is significantly larger than the others, 
we divided it by a factor of 10 in order to show it in the 
same plot. 

Models PI, P3 and P4 constitute a crude numerical 
resolution study for these models. In order to character- 
ize the effect of numerical resolution on our results we 
choose the run PI with numerical resolution of 100'^ as 
a baseline model. With tangled azimuthal field struc- 
ture and characteristic coherence lengths of ri = 0.1 and 
r2 = 0.087 this run is a good choice for resolution study, 
because the resolution effects are potentially more severe 
than in the runs with large scale, loosely tangled field. 
We carried out calculations equivalent to PI, at higher 
resolution of 200'^ (run P3) and lower resolution of 64'^ 
(run P4). Interestingly, we find that stepping up the nu- 
merical resolution leads to a non-monotonic evolution of 
the collapse timescale (Fig. 16). This behavior can be 



explained by the effect of resol utio n on magnetic field 
orientation (right panel of Fig. 16 1. Higher resolution 



captures the finer, radial structure of tangled field lines 
in addition to the dominant azimuthal component, re- 
sulting in the lower initial value of {Ob)- This results 
in a highest rate of heat conduction in model P3 (0.04, 
0.054, and 0.01 of the Spitzer value for PI, P3, and P4, 
respectively) and a more rapid reorientation of the field 
lines due to HBI when compared to moderate-resolution 
model PI. The consequence of rapid HBI evolution in 
model core P3 is that it also collapses before the core in 
PI. On the other hand, heat conduction in model P4 
is initially at such a low level that its evolution proceeds 
fairly similar to a pure radiative collapse and thus, occurs 
on a shorter time scale that in models PI and P3. Pro- 
gressing from the simulation with 64"^ resolution towards 
the higher end, the discrepancy in overall evolution time 
of the core decreases from ^ 2 Gyr to ~ 1 Gyr, indi- 
cating convergence. Based on this we estimate that core 
evolution in most affected models in our calculations can 
be prolonged by ~ 1 Gyr. 

4. DISCUSSION 



Heat conduction is in principle capable of providing 
a large fraction of the energy necessary to prevent the 
core collapse in clusters of galaxies — indeed, the po- 
tential importance of conduction has been discussed ever 
since the realization that thermal collapse was a prob- 
lem for ICM atmospheres (Binney & Cowie 1981). How- 
ever, attempts to solve the cooling flow problem with 
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Fig. 14. — Snapshots from PI model of the Perseus cluster showing the distribution of velocity magnitude at 500 Myr (left) and 3.5 Gyr 
(right) after the beginning of the simulation. Figures show the slice through the xy-plane of a cluster center. A velocity unit corresponds to 
820kms~^, implying that highest velocities are of the order of 30kms^^. The apparent symmetry in velocity distribution is a consequence 
of initial symmetry in magnetic field geometry. 
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previous studies have highlighted fine-tunin g problems 
( [Nulsen et al |1982| [Bregman fc D avid 1988). Under the 
assumption tliat a tangled magnetic held essentially acts 
as a scalar suppression factor for isotropic conduction, 
and that the suppress ion factor / is a fixed parameter, 



2 4 6 8 10 

Time (Gyr) 

Fig. 15. — Evolution of luminosity components in model PI 
measured within a sphere of radius of 20 kpc. Each component is 
marked with a different line: anisotropic heat conduction (solid), 
radiative (dashed), convective (dash-dot), advective (dash-Sdot), 
and a fiducial value of unbridled heat luminosity at the Spitzer 
value (thick, long-dashed). The heat flux at the full Spitzer value 
was divided by a factor of 10 in order to show it alongside of other 
components. Negative fiuxes are inflowing and vice versa. 



thermal conduction face several challenges. Because of 
a steep temperature dependence of Spitzer conductiv- 
ity, heat conduction is expected to play a lesser role in 
low mass clusters and gro ups o f galaxies, where most 
of the gas is below 5keV ( Best et al.^2007; Zakamska, 
fc Narayan 2003). Even unbridled Spitzer thermal con- 
duction cannot offset radiative loses in these low-mass 
systems. Hence, heat conduction cannot be a single so- 
lution to the cooling flow problem across the full range 
of masses, and additional sources of heat are needed to 
stabilize cooling flows in at least these lower mass sys- 
tems. Even in hot systems where conduction is potent. 



Kim fc Narayan ( 2003 1 show that / must be flne-tuned in 
order to obtain a thermal quasi-equilibrium resembling 
real cooling core clusters. Even then, this equilibrium is 
unstable with a growth time ~ 2 — 5 Gyr. 

Our results indicate that presence of the HBI even fur- 
ther undermines the role of heat conduction because of its 
tendency to stifle conduction by rearranging the lines of 
magnetic fleld to be orthogonal to the temperature gra- 
dient. We find that only systems with a large disparity 
in characteristic time scales, tcooi/icond ^ 10^, can avoid 
the cooling catastrophe by evolving to isothermal state. 
Systems reminiscent of real clusters, with icooi/^cond ^ 25 
evolve towards thermal collapse on time scales < lOicooi- 
For many cooling core clusters (such as Perseus) this is 
significantly shorter than a Hubble time, hence render- 
ing heat conduction alone incapable of rescuing the core 
from collapse. 

That said, we do find that conduction can significantly 
increase the time to thermal collapse, i.e., a significant 
energy is conducted into the core before the field-line 
re-orientation by the HBI seals it off and enables it 
to collapse. At this point, we must acknowledge the 
missing ingredients in our idealized cluster models. A 
quiescent cluster core can be significantly disturbed by 
cosmologically-induced dynamics (e.g. sub-cluster merg- 
ers) and/or powerful outbursts from a central AGN. Ei- 
ther of these events will disrupt the azimuthal nature of 
the field and may partially reset the evolution being de- 
scribed by our models. Thus, while thermal conduction 
seems unable to stabilize cooling cores alone, it may still 
be the primary agent heating cluster cores provided that 
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Fig. 16. — Effect of numerical resolution on evolution of the core temperature (left) and (Sg) (right) in PI (solid), P3 (dashed), and P4 
(dash-dot) runs carried out at resolution of 100"^, 200'^, and 64^, respectively. 



it is enabled by subcluster mergers or AGN. Note that 
this is a very different role than that normally envisaged 
for AGN in cooling cluster clusters; if correct, AGN are 
stirrers and not heaters. 

An interesting property of our models is the absence 
of strong turbulent motions. With plasma motions at ^ 
10—30 kms~^ modeled collapsing cores are kinematically 
quiescent in comparison with realistic cores, which ex- 
hibit turbulent velocities in the range 100 — 200kms~^ 



ove r comparable spatial scales (Conselice et al. 2001 
[Hatch et al.||2006| . This highlights the necessary role 
of the central AGiSf and cluster mergers as stirrers of the 
intercluster plasma. 

Both observational and theoretical studies of the cool- 
ing problem have established that the mechanism respon- 
sible for the stability of cooling cores has to be gentle, 
spatiall y distributed, and physically reminiscent of dif- 
fusion ( |Kim fc Narayan||2003| [Fabian et al.||2003j ). This 
is an intrinsic property of heat conduction but is more 
difficult to understand in pure AGN heating scenarios. 
For example, simple shock heating by AGN jets fails to 
stabilize cooling flows in pur e hydrodynamic simulations 
( Vernaleo fc Reynolds[[2006 ) unless a high degree of tur - 
bulent heat diffusion is postulat ed ( Briiggen et al.|2009 l. 
In addition, Nagai et al. (20071 find that the entropy of 
the gas in cooling cores is surprisingly low - at 10s of 
keV/cm^ the gas is on a brink of catastrophic collapse 
and yet maintained in that state over a large fraction 
of the Hubble time. This presents a challenge for any 
mechanism which would heat the gas via strong turbu- 
lent motions or shocks, as it is expected to result in a 
much higher value of entropy than that observed. The 
implication of this discovery is that an unusually deli- 
cately balanced feedback mechanism operates in cooling 
cores and that heat conduction in combination with ad- 
ditional mechanisms may play an important role. 

Another implication of the ubiquity of HBI is that clus- 
ters with cool cores may exhibit a characteristic magnetic 
field structure, where field lines are wrapped around the 
core and have predominantly azimuthal structure. Al- 
though our model clusters never completely forget about 
their initial magnetic field (i.e., T-models and A- models 
can be recognized even at late time), the wrapping of field 
lines into spherical shells is ubiquitous. The instability 



drives the field lines to almost orthogonal position with 
respect to the temperature gradient, {9b) ~ 80 — 90°. 
This has several consequences. Firstly, it creates ex- 
actly the kind of magnetic field geometry envisaged in 
models of AGN blown bubbles that invoke "magnetic 
draping" (Ruszkowski et al. 2007). In these models, a 
bubble of relativistic plasma which is buoyantly rising in 
the ICM atmosphere is stabilized to Rayleigh- Taylor and 
Kelvin-Helmholtz instabilities by a layer of strong mag- 
netic field that has been swept up on its leading edge. It 
is suggested that this mechanism is responsible for the 
surprising stability of "ghost cavities" in the ICM of the 
Perseus cluster, and is an alternative to models that in- 
voke plasma viscosity (Reynolds et al. 2005) . Secondly, 
this characteristic field line geometry leads to the most 
direct observable predictions/manifestations of the HBI 
in the ICM. Magnetic fields in the ambient ICM can be 
probed via the Rotation Measure (RM) of polarized radio 
sources in the background. 

The RM is an integral measurement of the effect that 
magnetic fields impart on the orientation of the plane of 
an electromagnetic wave traveling through the plasma. 
Hence, an azimuthal field structure would produce a 
characteristic radial dependence in the observed RM: it 
would vanish at the projected cluster center, and peak 
at some projected radius within the cooling core. We 
calculate this effect for our Perseus-like models using the 
expression RM= 812radm^^ J n^'B ■ dl, where rif. is the 
electron number density in units of cm~'^, B is the vec- 
tor of magnetic field in units of fiG, and 1 is the depth 
of the magnetic screen in kiloparsecs. Figure [T7| shows 
an illustrative calculation of the RM for model-P2. We 
find that in Perseus-like models the magnetic field with 
mean strength of (B) ^ 0.1/iG produces RM of order 
~ lOOradm"^. The HBI-wrapped field also produces a 
characteristic spatial structure to the RM map. Future 
observatories with the radio instruments such as Square 
Kilometer Array should provide the sensitivity needed 
to identify and study background sources with sufficient 
areal density to map out these RM patterns. Note how- 
ever that our calculation accounts for the RM arising 
from the core of a cluster but does not take into account 
the effect of magnetic structures and plasma in the outer 




Fig. 17. — Rotation measure maps calculated for model P2 for a line of sight along x-axis. The colorbar indicates the RM intensity which 
ranges from (—max, max), where max= 150, 110,95, 120, 170,265 radm^^^ foj, g^ch panel in the time sequence, respectively. 



parts of a cluster. Therefore, if the cluster field strength 
is dominated by the the magnetic fields in the core, this 
effect may be measurable. 

5. CONCLUSIONS 

We have performed simulations of cooling cores in 
clusters of galaxies with the aim to study the role of 
anisotropic heat conduction and recently discovered HBI 
on the thermodynamic evolution of cores. Our models 
focus on the base state of cluster cores and do not take 
into account physics of the central AGN or dark matter 
or ICM substructure resulting from hierarchical merg- 
ing of subcluster units. We explore the parameter space 
of cooling and conduction timescales as well as different 
magnetic field configurations. We summarize our most 
important results here: 

• The parameter space of modeled systems can 
be divided into three distinct groups of mod- 
els where collapsing cores occupy the parameter 
range icooi/^cond ^ 25 and isothermal cores exhibit 
icooi/icond ^ 10^. The range between these two 
groups is occupied by a class of systems which fi- 
nal state is determined by the configuration of mag- 
netic field lines. 

• The efficiency of heat conduction along the lines of 
magnetic field ranges between ~ 10~^ — 0.2 of the 



Spitzer value, depending on the initial structure of 
magnetic field and evolutionary stage of the core. 

• Modeled conducting cores that correspond to real 
clusters, exhibit evolution towards core collapse 
slower by a factor of '--^ 2 — 10 with respect to 
the time scale for a pure radiative collapse. For 
many cooling core clusters (such as Perseus) this 
is significantly shorter than a Hubble time, imply- 
ing that heat conduction alone cannot rescue the 
core from collapse. The extent to which core col- 
lapse is postponed in our models is a function of 
the initial magnetic field topology, where systems 
with higher values of tcooi/^cond ratio and field con- 
figurations amenable to conduction result in longer 
collapse time scales. 

• Magnetic field lines in cores in which HBI is ac- 
tively operating rearrange themselves in such way 
that final value of {0b) ~ 80° - 90°. In con- 
trast, cores that are not characterized by HBI have 
{9b) ~ 60° — 70°. If seen in observations, the align- 
ment of field lines across the global temperature 
gradient would provide a strong evidence for exis- 
tence of MHD instabilities in clusters, which have 
so far only been considered theoretically. 

• Heat flux buoyancy instability operating in mod- 
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eled cores is not characterized by a strong turbu- 
lence and results in kinematically quiescent cores. 
A related consequence of this property is that MHD 
turbulence driven by the heat flux instability does 
not drive a convective thermal flux suflicient to reg- 
ulate the core collapse. This finding is interesting 
in light of observational and theoretical studies of 
the cooling problem which have established that 
the mechanism responsible for the stability of cool- 
ing cores has to be gentle and spatially distributed 
in order to agree with the low levels of entropy ob- 
served in some cooling cores and implies that heat 
conduction in combination with additional mecha- 
nisms may play an important role for the cooling 
problem. 

The prospects for future work are numerous and they 
stem from the need to design more realistic models of 
cooling cores and also to understand the interplay of 
MHD instabilities with other mechanisms operating in 
clusters. Future studies will be needed to address the 
evolution of cool core clusters within the cosmological 
framework and consider the effects of minor and major 
mergers on gas mixing, kinematically and MHD driven 
turbulence, field line tangling, AGN fueling, and the 



overall thermodynamic state of a cluster. On the other 
hand they will also need to disentangle the effects of 
plasma transport processes and magnetic fields via sim- 
ulations that incorporate full MHD as well as viscosity, 
thermal conduction, and cosmic ra y diffusion. 



Note that in parallel to this work Parrish et al. ( 2009 1 



carried ou t an independent st udy of cool cluster cores 



with HBI. Parrish et al. ( 2009 1 focus on a portion of the 
parameter space occupied by real clusters and consider 
the role of central entropy and AGN heating for the final 
thermodynamical state of a cluster core. Given different 
approaches, the two studies complement each other and 
in the areas of overlap arrive to similar conclusions. 



We thank Sean O'Neill and Eve Ostriker for useful and 
insightful discussions. Numerical simulations of MHD 
instabilities in clusters were carried out on Deepthought, 
the University of Maryland high performance comput- 
ing cluster and on jxb at the Ecole Normale Superieure. 
TB thanks the UMCP-Astronomy Center for Theory 
and Computation Prize Fellowship program for support. 
CSR acknowledges support from the Chandra Theory 
and Modeling Program under grant TM7-8009X. 



REFERENCES 



Asai, N., Fukuda, N., & Matsumoto, R. 2004, ApJ, 606, L105 

Balbus, S. A. 2000, ApJ, 534, 420 

Balbus, S. A., & Reynolds, C. S. 2008, ApJ, 681, L65 

Best, P. N., von der Linden, A., KaufTmann, G., Heckman, T. M., 

& Kaiser, C. R. 2007, MNRAS, 379, 894 
Binney, J., & Cowie, L. L. 1981, ApJ, 247, 464 
Bregman, J. N., & David, L. P. 1988, ApJ, 326, 639 
Briiggen, M., Scannapieco, E., & Heinz, S. 2009, MNRAS, 395, 

2210 

Conselice, C. J., Gallagher, J. S., Ill, & Wysc, R. F. G. 2001, AJ, 
122, 2281 

Ettori, S., & Fabian, A. C. 2000, MNRAS, 317, L57 

Fabian, A. C., Sanders, J. S., Taylor, G. B., Allen, S. W., Crawford, 

C. S., Johnstone, R. M., & Iwasawa, K. 2006, MNRAS, 366, 417 
Fabian, A. C., Sanders, J. S., Crawford, C. S., Conselice, C. J., 

Gallagher, J. S., & Wyse, R. F. G. 2003, MNRAS, 344, L48 
Fabian, A. C, Voigt, L. M., & Morris, R. G. 2002, MNRAS, 335, 

L71 

Hatch, N. A., Crawford, C. S., Johnstone, R. M., & Fabian, A. C. 

2006, MNRAS, 367, 433 
Kim, W.-T., & Narayan, R. 2003, ApJ, 596, 889 
Nagai, D., Kravtsov, A. V., & Vikhhnin, A. 2007, ApJ, 668, 1 
Narayan, R., & Medvedev, M. V. 2001, ApJ, 562, L129 



Nulsen, P. E. J., Stewart, G. C, Fabian, A. C, Mushotzky, R. F., 
Holt, S. S., Ku, W. H.-M., &c Malin, D. F. 1982, MNRAS, 199, 
1089 

Parrish, I. J., Quataert, E., & Sharma, P. 2009, submitted to ApJ 
Parrish, I. J., Stone, J. M., & Lemaster, N. 2008, ApJ, 688, 905 
Parrish, I. J., & Quataert, E. 2008, ApJ, 677, L9 
Parrish, I. J., & Stone, J. M. 2007, ApJ, 664, 135 
Parrish, I. J., & Stone, J. M. 2005, ApJ, 633, 334 
Quataert, E. 2008, ApJ, 673, 758 

Sharma, P., & Hammett, G. W. 2007, Journal of Computational 

Physics, 227, 123 
Spitzer, L. 1962, Physics of Fully Ionized Gases, New York: 

Interscience (2nd edition), 1962 
Stone, J. M., Gardiner, T. A., Teuben, P., Hawley, J. F., & Simon, 

J. B. 2008, ApJS, 178, 137 
Tozzi, P., &; Norman, C. 2001, ApJ, 546, 63 
Vernaleo, J. C, &; Reynolds, C. S. 2006, ApJ, 645, 83 
Vikhlinin, A., Markevitch, M., & Murray, S. S. 2001, ApJ, 551, 160 
Voigt, L. M., Schmidt, R. W., Fabian, A. C, Allen, S. W., & 

Johnstone, R. M. 2002, MNRAS, 335, L7 
Zakamska, N. L., & Narayan, R. 2003, ApJ, 582, 162 



